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ABSTRACT 

In gravitationally stratified fluids, length scales are normafly much greater in the horizontal 
direction than in the vertical one. When modelling these fluids it can be advantageous to use 
the hydrostatic approximation, which filters out vertically propagating sound waves and thus 
allows a greater timestep. We briefly review this approximation, which is commonplace in 
atmospheric physics, and compare it to other approximations used in astrophysics such as 
Boussinesq and anelastic, finding that it should be the best approximation to use in context 
such as radiative stellar zones, compact objects, stellar or planetary atmospheres and other 
contexts. We describe a finite-difference numerical scheme which uses this approximation, 
which includes magnetic fields. 
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1 INTRODUCTION 

In magnetohydrodynamical (MHD) simulations, a set of partial dif- 
ferential equations is numerically integrated forwards in time. This 
is done in stages, with the time increasing in small increments 
called the "timestep". When the basic MHD equations are used, 
the timestep is subject to a range of limits to do with the speed of 
propagation of information; all numerical schemes become unsta- 
ble if information is allowed to propagate further than some fraction 
of a grid spacing in one timestep. For instance, a simple Cartesian 
hydrodynamical code might have the following timestep: 



As well as the basic constant-density incompressible approxima- 
tion, in which sound and buoya ncy waves are both absent, the 
anelastic l iOgura & Phillipsl 1962j) and Boussinesq approximations 



At = min 
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(1) 



where Aa;, Ay and Az are the grid spacings in the three dimen- 
sions, u is the gas velocity, Cs is the sound speed and C is some 
dimensionless constant whose value will depend on properties of 
the numerical discretisation scheme, but might be for instance 0.5. 
This would ensure that no information can propagate more than 
0.5 grid spacings in any direction during one timestep. Other codes 
will have additional, similar restrictions from the propagation of 
Alfven waves and from diffusion; for instance, the sound speed Cs 
in the expression above might be replaced by the fast magnetosonic 
speed. 

By studying the context in which we wish to use simulations, it 
is often possible to make approximations in order to remove some 
modes of propagation of information, allowing a larger timestep. 



E-mail: jonathan@astro.uni-bonn.de 



are commonplace (see e.g. Lillvll 19961 . for a review and compari- 
son). Both are widely used in astrophysical hydrodynamics, for in- 
stance in studies of convection in planetary a nd stellar interiors (e.g. 
lBrowning|l200i ; IChen & Glatzmaiedl2005l) . They can be used in 
situations where the thermodynamic variables depart only slightly 
from a hydrostatically balanced background state, so for instance 
the density perturbation Sp/po <C 1. There are some other require- 
ments, such as that the frequency of the motions is much less than 
the frequency of sound waves and that the vertical to horizontal 
length scale ratio or the motion is not too large. The two approxi- 
mations are rather similar, the difference being that the Boussinesq 
approximation is used where the vertical scale of the motions is 
much less than the density scale height and where the motion is 
dominated by buoyancy. It can be shown that the continuity equa- 
tion, whose standard form is dp/dt + ^ ■ pu — 0, reduces to the 
forms V • pou = and V ■ u = in the anelastic and Boussi- 
nesq approximations respectively. The result of this is that sound 
waves are filtered out and the timestep is no longer restricted by 
their propagation. Buoyancy waves are still allowed. 

In atmospheric physics, it is common to make the approxima- 
tion of hydrostatic equilibrium (equation |8), in which we assume 
perfect vertical force balance. This is applicable in contexts where 
a constant gravitational field causes strong stratification, where the 
length scales in the vertical direction are much smaller than in the 
horizontal, and where the fluid adjusts to vertical force balance on 
a timescale much shorter than any other timescale of interest - on 
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the timescale of so und waves propagating in ttie vertical direction 
( lRichardsonlll922h . The consequence is that vertically propagat- 
ing sound waves are filtered out, as well as high frequency inter- 
nal gravity waves, and the z-component of the timestep restriction 
(equation [T]l can be removed. This is an obvious advantage in any 
situation where the vertical length scales present in the system are 
much smaller than the horizontal length scales and adequate mod- 
elling therefore requires that Az <^ Ax, Ay. 

The hydrostatic approximation reduces the number of inde- 
pendent variables in the system. For instance, in a system where 
the gas has two thermodynamic degrees of freedom, in 'raw' hy- 
drodynamic equations there are these two thermodynamic variables 
plus the three components of velocity. In the equivalent hydrostatic 
system the vertical component of the velocity is no longer inde- 
pendent, but calculated by integration of equation ((S); furthermore 
one of the thermodynamic variables is lost. In the astrophysical 
context we often want to model conducting fluids with magnetic 
fields. Note that although the hydrostatic approximation filters out 
vertically-propagating sound waves, vertically-propagating mag- 
netic waves are not entirely filtered; a magnetohydrostatic scheme 
is therefore of use only in the case where the plasma /3 is high, i.e. 
where the Alfven speed is much less than the sound speed. 

There are various ways in which the hydrostatic approxima- 
tion can be implemented, resulting in different sets of equations 
and independent variables. The vertical coordinate can be physical 
height, pressure, entropy or some c ombination of those (see e.g. 
(Kasahara 1974: Ixonor & Arakawj[l997. , for a review of coordi- 
nate systems). It turns out, for instance, that a change of the vertical 
coordinate from height z (as used in other systems) to pressure P 
simplifies the equations. This can be seen by noting that in hydro- 
static equilibrium, the pressure at any point is simply equal to the 
weight of the column of gas above that point and that each grid box 
(which has a constant pressure difference AP from top to bottom) 
will contain constant mass. The continuity equation therefore be- 
comes dux/dx + duy/dy + duj/dP — 0, where to = DP/Dt 
the full Lagrangian derivative, which has the same form as the fa- 
miliar incompressible equation V ■ u = where vertical velocity 
Uz = Dz/Df has been replaced by u. Entropy coordinates (also 
known as isentropic coordinates) are also commonplace, the main 
advantage being that the vertical 'velocity' T)s/Dt is small, a func- 
tion only of heating and cooling, which reduces numerical diffusion 
in the vertical direction. 

The best choice of vertical coordinate often depends on the 
desired upper and lower boundary conditions. In weather forecast- 
ing, for instance, it is necessary to have the lower boundary fixed in 
space. Using pressure coordinates , implementation of this i s chal- 
lenging. It is for this reason that iKasahara & Washington! ( Il967h 
produced a hydrostatic numerical scheme using height coordinates, 
but owing to advances in hybrid coordinate systems which allowed 
also for topographical features - mountain ranges and so on - this 
scheme never became popular. However, when magnetic fields are 
added, height coordinates z will be simpler than either pressure or 
entropy coordinates and may regain an advantage in some contexts. 

In more astrophysical contexts, such as neutron star, stellar 
or planetary atmospheres, we may want a lower boundary fixed in 
space, and to be more precise, fixed at a particular height (unlike 
in the terrestrial context, mountain ranges and so on need not be 
included). It is often desirable to have the upper boundary fixed in 
pressure, if the temperature, and therefore also the pressure scale 
height, varies by a large factor. For instance, during X-ray bursts 
on neutron stars the temperature increases by about a factor of ten 
so that an upper boundary fixed in space would mean insufficient 



resolution of the relevant layers in cold areas, and extremely low 
densities and high Alfven speedsQ Entropy coordinates are unsuit- 
able since convection may appear, and in any case the entropy of 
a co-moving fluid element is expected to change rapidly, removing 
any advantages of this system. In this context, therefore, the natu- 
ral choice is the a-coordinate sy stem, a pressure-related coordinate 
first proposed bv |PhilliDsl ( ll957l) . 

In section[2]we present the basic equations, before describing 
the finite-difference numerical method in more detail in section [3] 
presenting simple test cases in section|4]and summarising in section 

m 



2 BASIC EQUATIONS 

In this section we describe the a-coordinate system and how mag- 
netic fields are incorporated. 

First of all, we describe the standard MHD equation^ and 
then go on to the additional equations coming from the hydrostatic 
approximation. Writing down the horizontal part of the velocity as 
u, the horizontal part of the momentum equation is 



= -VhP + ;^[V X B X Blh + Ff 

Dt A-K 



(2) 



where the subscript h denotes the horizontal component of a vector, 
and D/Dt = d/dt + u ■ Vh + itzC^/Sz is the Lagrangian deriva- 
tive. F^'^'^ is the viscous force per unit volume. Writing down the 
continuity, energy and induction equations and the equation of state 
(i.e. the perfect gas law), we have 



CP 



m = -Vh-(pu)-^(pn. 

DT _ IDP 

Dt ^ pDt 

^ = Vx(uxB - cE"^'^ 
dt ^ 

P = pRT, 



(3) 
(4) 

(5) 
(6) 



where Uz is the vertical component of the velocity, Q is the heating 
rate per unit mass (including that from heat conduction), cp is the 
specific heat at constant pressure, R is the gas constant (the univer- 
sal gas constant divided by the mean molecular weight of the gas 
in question) and other quantities have their usual meanings. In this 
system the vertical coordinate is a quantity a related to the pressure 
P and the pressure at the top and bottom boundaries Pt and Pb by 

_ P - Pt 



where 



P. 



(V) 



One of the important main features of this scheme is that here the 
upper boundary is fixed at pressure Pp (where cr = 0) rather than 
being fixed at a particular height. The lower boundary at cr = 1 
is fixed in space. Now, the standard set of MHD equations would 
also contain the z-component of the momentum equation (|2j but in 
the hydrostatic approximation we instead assume that the pressure 
gradient and gravity are always in perfect balance, i.e. we have the 
equation of hydrostatic equilibrium 



dP 
dz 



= -9P, 



(8) 



^ It is for this reason that Boussinesq and anelastic schemes are unsuitable 
here, since they cope with only small variations about a constant reference 



state. 

2 



We use c.g.s. units throughout. 
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neglecting the vertical component of the Lorentz force - which is 
much smaller than the pressure and gravity forces in this high-/? 
regime. 

In this scheme the fundamental variables which are evolved in 
time are the horizontal part of the velocity u, the temperature T and 
the pressure difference P*, plus optional quantities such as heating 
rate Q which we include here. Note that while T{x, y, a) and other 
variables are three-dimensional, P, (x, y) is only two-dimensional. 

The density, or rather its inverse a = 1/p, is calculated from 
the equation of state (EOS), which in the ideal gas case is 

RT 



(9) 



This can easily be modified to more complex EOS, for example 
around the transition between ideal gas and degenerate (electron) 
gas. As before, the vertical component of the momentum equation 
is replaced by the equation of hydrostatic equilibrium (|8) which 
takes the form 



dz 

gr — = -aP., 
oa 



(10) 



which we integrate from the lower boundary upwards to give height 
z and potential (p: 



— gz — Pt a da' 



(11) 



In some sense the quantity P, can be considered a 'pseudodensity' 
as it takes the role of density in the equations - the continuity equa- 
tion is: 



(12) 



where Vo- represents the gradient at constant a (unlike Vh which 
is at constant height). This equation has an obvious similarity to 
the familiar form dp/dt + V • (pu) = 0. The equivalent of the 
vertical component of the velocity is cr = Da/Dt. Note that P, is 
not a function of o and so comes outside of the derivative in the 
third term above. Given the boundary conditions that ct = at both 
upper and lower boundaries, we can integrate equation (1121 from 
cr = downwards to give 

dP, 



a— h/ + P.o- = where /; 

ot 



V„-(P*u)do-', (13) 



so that integrating to the lower boundary a — 1 gives a predictive 
equation for P, : 

dP, 

-or = '''' 

Substituting this back into equation il3\ allows calculation of the 
vertical velocity a: 



(15) 



Taking the Lagrangian derivative of the definition of a (equa- 
tion|7]l and multiplying by P, gives 



„ . DP DP. 

Pff^ ~ — , 

Dt Dt ' 

where the Lagrangian derivative is defined in the usual way 

Da „ . d 

\jt at oa 
which we can also use to express the time derivative of P« as 

DP. dP, 



Dt dt 



+ u<,=i ■ VP. 



(16) 



(17) 



(18) 



Substituting this, equations il4l and dlSt into equation ( 1161 ) gives 
us an expression for DP/Dt: 



DP r 
— = au....VP.-/ 



(19) 



which we need to evaluate the time derivative of temperature from 
the thermodynamic equation 



DT DP 
Cp-^ — = Ce— — 

Dt Dt 



(20) 



where Q, the rate of heating per unit mass, includes all heating, 
cooling and conductive terms. 

What remains now is the (horizontal part of the) momentum 
equation: 



Du 

Dt 



era VP, - 2n X u + Fvisc 



(21) 



It is not immediately obvious how best to go about adding 
magnetic fields to this scheme, since calculating real-space gradi- 
ents in the vertical direction necessitates first calculating the gra- 
dient 8/ da and then multiplying by da/dz = —g/aP,, and also 
because the grid points themselves are moving in the vertical direc- 
tion. 

The best way is to start by making a switch of the independent 
variable from B — {Bx, By B^) to 



B* 



dz „ dz 



B^—^B: 



(22) 



the equations can be significantly simplified, mainly since the time 
derivative {dz/dt)„ is not required; what we do require is just the 
derivative dz /da, which we have already from equation dlOt . Fur- 
ther defining V* = (d/dx,d/dy,d/da) and u* — {ux,Uy,a), 
we have 



5B' \ , , , . 

— ) =V x(u xB 



e: 



(23) 



This system ensures conservation of flux. Details of the viscous part 
of the electric field are given in section |3^.4.3l 

The Lorentz force is calculated by first dividing the x and y 
components of B* by dz/da to find the actual magnetic field B, 
then finding the current J from V x B in the usual way (where 
the vertical derivatives are of the form (da/dz)d/da, and simply 
taking the cross-product of the current with the magnetic field. 



3 NUMERICAL IMPLEMENTATION 

We now describe how the above equations are integrated numeri- 
cally, describing the grid, timestepping and the diffusion scheme. 



3.1 Numerical grid 

The grid is staggered, which improves the conservation proper- 
ties of the code and is worth the modest extra computational ex- 
pense; different variables are defined in different positions in the 
grid boxes. The horizontal components of the velocity are face- 
centred, defined half of one grid spacing from the centre of each 
grid box (see Fig.|2j, whilst all other variables are defined either in 
the centre of the grid box or vertically above/below it. Temperature 
and any other thermodynamic variables such as Q, and potential <j!> 
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are defined at the centre of each grid bojjfl P* is a function of x and 
y only and is defined in the centre of the grid box, not displaced in 
the X 01 y directions. Bx, By and Bz are face-centred, defined half 
of a grid spacing from the box centre in the x, y and z directions 
respectively. 

At various times while evaluating the time derivatives in the 
partial differential equations given above, it is necessary to find the 
spatial derivatives of various quantities, to evaluate a quantity at a 
position other than that where it is defined (e.g. half a grid spac- 
ing displaced in some direction) and to integrate a quantity in the 
vertical direction. The derivatives, interpolations and integrals are 
evaluated to fifth, sixth and fifth order respectively, meaning that 
the values of the given quantity at six grid points are used to calcu- 
late the re quired quan tity /derivative/integral at each grid point (for 
details see Eei3ll992h . 

For instance, if a quantity / is defined displaced half a grid- 
spacing in the x-direction from the centre of the grid box and its 
value is required at the grid-box centre the value is calculated thus: 



75 



256 ^'^'+5 



(24) 



and if the spatial derivative with respect to x is required at the same 
location, it is calculated thus: 



A f' 225 
Ax/, = —(U 



- — (f 



(25) 



Note that this method of calculating interpolations and derivatives 
is also used in the 'stagger code' llNordlund & GalsgaardI [l995l : 
iGudi ksen & Nordlundl llOOsI) . 

In addition to this, however, we need to integrate various quan- 
tities in the vertical direction - for instance to integrate a body- 
centred quantity / (ind ices 1/2, 3/2, etc.) from the upper boundary 
(where coordinate a = 0) downwards to an arbitrary position a, 
and to return the result at face-centred locations (indices 0, 1, 2 
etc.) we first perform a first-order integration 



rlst 



(26) 



and then increase the order with the following operation to add 
parts onto either end of the integration: 



= A 



'''+b{f_.-f.) + c(f_. 



) + (27) 
0.00484768 Act and d = 



where b = 0.0543134Acr, c = 
0.000379257Acr. 

In addition, a body-centred integrand must sometimes be inte- 
grated from the upper boundary (ct = 0) to body-centred positions 
As before, a first-order result is obtained first: 



rlst r-ist I Act 



. , , w,-^+,/;+i) (28) 

and a higher order result for just the point half a grid-spacing from 





cdz3 Bz, li = ^^..-^ 




^ T. Q. z. 
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i+l/2 




, T. Q, z. 
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N- 1/2 



^Bz,.. 



Figure 1. The grid, showing the positions within each grid box at which the 
fundamental quantities are defined: T, centres of grid boxes, dark circle; 
Vx, o" and y centred, x face centred, open circtes; Vy, a and x centred, y 
face centred, open circles. P* (x, y) is defined at the same x and y locations 
as T; any additional thermodynamic vaiiables such as the heating rate Q, 
are defined in the same locations as T. Bx is defined at the same positions 
as Vx ; By at the same as Vy . Bz and a are x and y centred, but face cen- 
tred in the cr direction. The an'ows represent positive directions for all the 
respective vector Kittie's, but &: it has the opposite sign. See text for defini- 
tions. Above and below the central box, the upper and lower boundaries are 
also shown. 
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other variables 



Figure 2. The so-called C-grid. Projection on the horizontal plane. Veloc- 
ities are face-centred (dark circles); the open circles in the centre of each 
grid box represent all of the other non-magnetic variables. Note that this 
diagram gives no information regarding the vertical positions. 



the upper boundary: 

II = tsi/i + bifs + ci/s 

2 2 2 2 

+a2/_i +&2f_3 +C2f_5. (29) 

2 2 2 

where the coefficients have the values 



^ The staggering in the vertical is known as the Lorenz grid iLorenzl l960h. 
The alternative is the Chamey-PhiUips grid iChamev &PhilliDslll953h 
where these quantities are face-centred, displaced half a grid spacing in the 
vertical from grid-box centi'e. 



ai = 0.41410590Acr aa 
bi = -0.036306424Acr 62 
ci = 0.0041015625Acr C2 



0.14283854Acr 
-0.028276910Acr 
0.0035373264Acr 
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These are used to produce the final result: 



(30) 



with a = (1/2)A(T, h = -(41/720)Acr and c = (ll/1440)Ao-. 
In other situations, it is convenient to perform this integration in the 
other direction, in which case the equivalent can be done, running 
the loops in the opposite order. 

It can be seen above that for interpolations, derivatives and in- 
tegrations, the values of quantities are required beyond the bound- 
aries of the computational domain; this is described below in sec- 
tion[331 



3.2 Boundaries 

In the horizontal directions, the simplest boundaries to implement 
are obviously periodic boundaries, but there is the possibility of 
switching one of the two directions to some other boundary condi- 
tion. Another possibility is to have 'mirror' conditions in one direc- 
tion, in order to avoid modelling the same thing twice in symmet- 
ric configurations. These conditions are symmetric in T and other 
scalars, and antisymmetric in the perpendicular component of the 
velocity. 

In the vertical direction, periodic boundaries are impossible, 
so symmetric conditions are used in T (and other thermodynamic 
variables) and the parallel components of velocity. These represent 
all of the independent variables, as the vertical component of ve- 
locity is a derived quantity. The zero condition (and antisymmetry) 
for the vertical velocity ensures that no mass can flow across the 
vertical boundaries. 



3.3 Timestepping 

The tim estepping uses a third-order Low-Storage Runge-Kutta 
scheme ( IWilUamsonlll98(]|) . which in practice means that during 
each timestep the time derivatives on the left-hand sides of the par- 
tial differential equations are evaluated three times, with three dif- 
ferent values of the quantities on the right-hand sides. If the quan- 
tities /, g, etc. are to be evolved in time from time t = to, at which 
the quantities have values /o, go, etc., each time step consists of the 
following steps (just including variable / for brevity): 

• Calculates time derivatives /q from /o 

• Finds timestep At according to various Courant conditions 
(see below) 

• Evaluates new values f\ 

h = /o + Atbi/^ 

• Evaluates new time derivatives f[ from /i and averages with 
previous one 

• Updates new values for second step 

/2=/l+At62/(.5 

• Calculates time derivatives /2 from /2 and averages with pre- 



vious ones 



• Calculates final time step values /s 

/3 =/2+Atfe3,f2.5 

and so /a is the value at the end of the timestep. The coeffi- 
cients are a2 = -0.641874, as = -1.31021, bi = 0.46173, 
&2 = 0.924087 and 63 = 0.390614. The advantage of this type 
of scheme is that the results from the previous evaluations need 
not be stored in memory, therefore making the code faster and less 
demanding in terms of memory usage. 

The time step is limited by the various Courant conditions 
given by the different quantities that are evolved. First defining 

(31) 

(32) 

(33) 

(34) 
(35) 



As = min(Aa;, Ay, Aa/ 

A = (6.2x3x0.75/2)^ 
As^ 

we have the following Courant parameters 
At 



max(cs + |u|) 



As 



At 



]_dT 
T~dt 

Amax[3 max(i^), max(!/s 



Cu is the limiting factor from velocities (both physical and sound 
velocity Cs), Cp the one from temperature changes, while Ci, is the 
one from kinetic diffusion (see sec. I3.4t . Finally, the timestep is 
fixed according to 

CAtAt 



At : 



(36) 



max(Cu, C„, Cp) 
where CAt is some numerical factor, which we generally set to 0.3. 

3.4 Diffusion 

In this section we describe how the scheme handles kinetic, thermal 
and magnetic diffusivities. 

The code includes both pure physical diffusion as well as a 
'hyperdiffusive' scheme, designed to damp structure close to the 
Nyquist spatial frequency while preserving well-resolved structure 
on larger length scales. Often it is possible, by assessment of their 
relative magnitudes (also treating the horizontal and vertical direc- 
tions separately) to switch off one or the other. The kinetic, thermal 
and magnetic diffusivities v, n and rj all have dimensions of length 
squared over time. 

It can be shown that with the physical 'textbook' diffusion 
equations, the computational demands sometimes become pro- 
hibitively expensive. For instance, using a kinetic diffusivity u high 
enough to handle shocks would result in unacceptable damping of 
low-amplitude sound waves unless one were able to use an un- 
realistically high resolution. Moreover, in constructing the fluid 
equations we made approximations based on the assumption that 
all relevant length scales in the system were much larger than the 
mean-free-path of the molecules, but the same fluid equations pro- 
duce shocks which physically have a thickness comparable to the 
mean-free-path, i.e. the equations have predicted their own invalid- 
ity. To allow shock handling without prohibitively high diffusivity 
in the entire volume, there is an additional viscosity in the vicinity 
of shocks. The two viscosities are one based on the sound speed 
and one on the divergence of the (horizontal) velocity - the latter is 
very negative in a shock. They are: 



V = !/i(cs + |u|) max(Aa::, Ay) 
Vs = VI smooth[max(— V ■ u, 0)] max(Aa::, Ay)^ 



(37) 
(38) 
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where both coefficients vi and i>2 are dimensionless; despite 1/2 
generally being much larger in value, the second term above is neg- 
ligible except in shocks. 

In addition to such a method of shock handling, many 
general-purpose codes use some kind of artificial diffusion scheme 
which can handle discontinuities and damp unwanted 'zig-zags'. 
Here, we use a 'hyper-dif fusive' scheme, based on that of 
iNordlund & GalsgaardI ( [l995h . where the diffusion coefficients are 
scaled by the ratio of the third and first spatial derivatives of the 
quantity in question, which has the effect of increasing the diffu- 
sivity seen by structures on small scales where the third derivative 
is high, damping any badly resolved structure near the Nyquist spa- 
tial frequency, while allowing a low effective diffusivity on larger 
scales. The way this works in practice is via diffusive flux opera- 
tors: 



where 
and 



d3. 



^y. ^ max(|d3,;+i|,|d3i|,|d3i-i|, ^^^^ 



max(|d/i+i|, |d/i 
U+i -/,_i)/Ax 
d,A+i-2d/, +d/,_i. 



|d/.-i 



(40) 
(41) 



These flux operators replace the derivatives of quantities on which 
the diffusion is operating, such as inside the brackets in equation 
i42\ below. For some kinds of diffusion we use these hyperdiffu- 
sive derivatives, and for other kinds we use the standard derivatives 
to give a more physical result. In addition, because of the very dif- 
ferent length scales and grid spacings, diffusion in the vertical and 
horizontal directions must often be treated differently. 



3.4.1 Kinetic diffusion 

Assuming that bulk viscosity is zero (a good approximation in 
monoatomic gases), the result of viscosity is to add the following 
viscous force (per unit mass) to the momentum equation (|2): 



1 d 

pdxj 



dui duj 



da 



+ 



dxi 



3^.. V 



(42) 



using Einstein summation notation. Furthermore, in many applica- 
tions it can be shown that the divergence term is very much smaller 
than the other terms, and we drop this term here, as it does not con- 
tribute in any case to numerical stability. Moreover, in many astro- 
physical applications including those for which this code has so far 
been used, the kinetic diffusivity is much smaller than the other two 
diffusivities and it is desirable to reduce the 'effective' viscosity as 
much as possible whilst preserving the stability of the code. For 
this reason, kinetic diffusion uses hyperdiffusive derivatives (de- 
scribed above) inside the brackets in equation ( 142) : the derivative 
outside the square brackets remains a standard high-order deriva- 
tive to preserve momentum conservation. However, the hyperdiffu- 
sive derivatives are used only with the standard viscosity (equation 
|37| l and the normal derivatives are used with the shock-viscosity 
(equation|38). 

The vertical and horizontal directions require different treat- 
ment. First, note that the vertical component of the viscous force 
is ignored, as we are not considering the vertical part of the mo- 
mentum equation. Second, all terms in equation J42) which contain 
the vertical velocity are dropped. Third, all derivatives with respect 
to the vertical coordinate must be scaled by a factor Aa/Ax or 



Aa/ Ay. The viscous force (per unit mass) is: 



1 

d_ 

dy 
d_ 
da 
j_ 



dx 
d 



d_ 

dx 

~2~ 



8*1 



dy 

P^vu„ 9* It, 



2 
d_ 

dy 
P,v 



dx 

d*Ux ^ d*u 
+ 



dux 
dx 



dx J ' 2 

P^Vs^a du. 



dux 
dy 



+ 



du^ 



P*u 



da 

d u- 



dy 



da 



(43) 



dy 



d*Uy d*Ux 



doL 



+ 



dy 



+ 



P*Vs ( dUy 

dx 



+ 



dux 
dy 



P*Wcr d*Uy 



PtUsUa dUy 



da 



da 



(44) 



where the asterisks with the derivative signify a hyperdiffu- 
sive differentiation as detailed above, and where is equal to 
Aa^ I min(Aa::, Aj/)^. These values of ensure not only that the 
units of the different terms are the same, but also that a given zig- 
zag structure is damped on the same number of timesteps, indepen- 
dently of the grid spacing in the three directions. Finally, note the 
role of P, , the pseudo-densit)0 - its presence in this way in the 
equations ensures conservation of momentum. 

In addition to the viscous stress, viscosity heats the fluid. The 
magnitude of this heating (per unit mass) is 



^^visc 



1 „ dui 



(45) 



where Sij is the viscous stress tensor, i.e. the contents of the square 
brackets in equations ( I43t to ( I44t . The index i is equal to x and y, 
but the index j is x, y and a. Units are taken care of by the presence 
of Vcr inside the Sxa and Sya parts of the stress tensor. 



3.4.2 Thermal diffusion 

The code includes two kinds of thermal diffusion: hyperdiffusive 
(similar to the momentum diffusion), and 'physical', both of which 
are simply added to the heating per unit mass Q. The hyperdiffusive 
thermal diffusion is 



?thcrm = — V ■ [P*Cp(KV*r 



(46) 



where the asterisk signifies a hyperdiffusive derivative as described 
above. The vertical derivatives are scaled with Aa'^ / Ax^ as before. 
The diffusivities k and Ks are calculated simply by multiplying v 
and i^s by a number, normally unity. 

It is sometimes desirable to have a larger thermal diffusion 
to model an actual physical process. To this end, the code also in- 
cludes a physical thermal diffusion, which is simply the same as 
equation i46t but with standard spatial derivatives. In the vertical 
direction, derivatives with respect to a are scaled with da /dz. In 
principle, there are also cross-terms originating from the fact that 
surfaces of constant a are not horizontal. Generally though, these 
terms are tiny and can be dropped. Also, the difference of length 
scales in the horizontal and vertical often means that the physical 
thermal diffusion in the horizontal direction is too small to provide 
numerical stability, and so a hyperdiffusive horizontal diffusion is 
required; likewise, in many applications the physical diffusion in 



* The right tenn should be P» /g, but g, being a constant, simplifies out of 
the equations. 
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the vertical direction is much larger than the minimum required for 
stability and the vertical hyperdiffusion can be switched off. 

3.4.3 Magnetic diffusion 

Finite conductivity gives rise to an extra electric field in the in- 
duction equation ([5} since in a medium with finite conductivity, an 
electric field in the co-moving frame is required to drive a cuiTent 
according to Ohm's law J = crE where here a is the electrical 
conductivity. Remembering that J — (c/47r)V x B and that the 
magnetic diffusivity is defined as r; = c^/(47ra), this extra field is 



-V X B. 



(47) 



The code contains a hyperdiffusive scheme rather like that de- 
scribed above for momentum diffusion. The electric current is cal- 
culated with the standard derivatives but the diffusivity rj is scaled. 
The expression for the electric field is 



E^'- = i^{r,h(J) + ,,.J} 



where the hyperdiffusive operator h is 



ha: = 7^ ( Ay^ 



9y2 



(48) 



(49) 



with corresponding values for the y and z components. 

The value of rj is given by multiplying v by some number, nor- 
mally unity. We determine rjs by a similar method to that used in 
determining Us, with the difference that we use the divergence not 
of the velocity field u but of the part of the velocity field perpen- 
dicular to the magnetic field, ux . 

The energy consequently lost from the electromagnetic field 
appears as heat, the so-called Joule heating given by (per unit mass) 



Qj 



= ij. 



(50) 



(51) 



3.4.4 Diffusion of other variables 

In addition to horizontal velocity, temperature and magnetic field, 
it is often necessary in the a-coordinate scheme to apply a hyper- 
diffusive scheme to the other main variable, P, . This works in ex- 
actly the same way as for the other variables, except that being 
two-dimensional it is somewhat simpler. The diffusivities v and Va 
are first averaged over a and multiplied by some number, normally 
unity, then a term is added to the partial differential equation I ll4t : 

dP* _ 

~dr ^ 

where the gradients are two-dimensional. This diffusion helps to 
damp unwanted zig-zag behaviour where it occurs. 

Any additional variables must generally also have some added 
diffusion: this works in exactly the same way as the hyperdiffusion 
on other variables. Normally the diffusivity can be lower, though, if 
the variable is a passive tracer with no feedback on other variables. 



3.5 Parallelisation and horizontal coordinates 

The code has been parallelised using OpenMP, which can be 
used on a shared-memoiy machine. Parallelisation for a distributed 
memory machine using MPI is planned in the medium term. 

Also planned for the medium term is an extension to spherical 
coordinates, with a view to modelling oceans and atmospheres on 
stars and planets. 



4 TEST CASES 

In this section we describe the numerical tests used to validate 
the code. We simulate the development of a shock from a wave, 
the Rossby adjustment problem and the Kelvin-Helmholtz and in- 
verse entropy gradient instabilities. We also test the propagation of 
Alfven waves and the Tayler instability for the magnetic field. 

All simulations have periodic horizontal boundary conditions. 
We assign the pressure via Pt, which is constant at all times, and 
P* (see equations 171 and [74l(. We also assign the temperature and 
the initial velocity fields. When initial prescriptions are better ex- 
pressed in terms of density, we assign temperature such that the 
right density is regained (see equation|9). 

One point has to be noted about the vertical coordinate in the 
figures. Given that we use the cr-coordinate system, neither phys- 
ical height nor gravity enters the equations (see equations II Hand 
I21t : only the product of the two is present. What is really relevant 
is the ratio between the scale height and the physical height of the 
model. What this means is that gravity g is essentially a free scal- 
ing factor which allows us to translate our simulations to different 
physical settings and the choices of gravity made here are arbitrary. 
However, when magnetic fields are added physical height becomes 
meaningful in its own right. 



4.1 Wave developing into a shock 

As a first test we start with a wave developing into a shock. These 
runs were performed in 2D, x and a. The experiment is set-up 
with a uniform temperature RT = 10^ erg g~^ and a perturba- 
tion in pressure (i.e. in P») such that the resulting perturbation in 
the height is sinusoidal 



5H/Ho = 0.1cos(27ra;/A) 



(52) 



where we used A equal to the extent of the domain (1 cm). 

The three resolutions used to check convergence were 50x50, 
100x100 and 200x200 (Fig. O. The wave develops into a shock 
after ~ 5 crossing times in all three simulations. We further use this 
set-up to check if the shock jump conditions are met (see below) 
and to do this we also run a new simulation of a strong shock at 
resolution 200x200 increasing the height perturbation amplitude to 
0.7 times the background value (Fig.|4j. 



4.1.1 Conserved quantities 

At this stage it is useful to check the conservation properties of the 
numerical scheme, in terms of mass, momentum and energy. 

We calculate the integrals as (remember that P./g is "equiva- 
lent" to p): 



Px / y ,tot 



E - 



P. If 

— dx dy da = — / P* da; dy (53) 

g g J 

P, 

Ux/y — da; dy da (54) 

ti? + uf. \ P, 

^ + CpT — da; dy da (55) 

^ J g 



In the last equation (see lKasaharal I974I . equation 5.18) (m^ +Uy ) /2 
represents the kinetic energy and cpT is the specific enthalpy. Note 
that there is no term for the gravitational potential energy - that is 
because as a whole, that energy is built into the enthalpy cpT. To 
see this, imagine 'inflating' the ocean from (close to) absolute zero 
whilst retaining the vertical ordering of each fluid element; for each 
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Figure 3. Top: Initial snapshot of the wave simulation 200x200 (section 
14. U . Units of the velocity color scale in cm s~^. Bottom: Snapshot at 
t = 2.44 X 10~^ s, first crossing time. Units of the velocity color scale 
in cm s~^. Supeiimposed is the initial height profile. The configuration is 
symmetrical with respect to the initial one. 



fluid element one needs just the eventual internal energy and the 
P dV work which is used in pushing the overlying fluid upwards, 
and the sum of these two is simply the enthalpy; the pressure of 
each fluid element remains constant during this process. 

We do not plot i\ftot conservation, because it is conserved to 
machine accuracy in all simulations. Fig. |5] shows the evolution in 
time of {E—Eo)/Eo for the three resolution simulations, where Eq 
is the initial energy. Energy is conserved to about one part in 10''. 
In these wave simulations the conservation of energy is essentially 
independent of resolution. Momentum conservation is looked at in 
section |43] since here the total momentum is zero and a fractional 
conservation is tricky. 








t = 4.34e-03 s 






Figure 4. The wave simulation 200x200, case of the strong shock (section 
14. It . Units of the velocity color scale in cm s"'^. The three frames are a 
time sequence at t = 0, 3.56 X 10"'^ and 4.34 X 10^^ s. Superimposed 
on the latter two frames is the initial height profile. The shock develops and 
is then reflected off the boundary. 



4.1.2 Wave speed 

We now study the velocity of the wave: considering only half the 
domain, since the horizontal velocity field Ux is zero at the bound- 
aries, before the nonlinear effects become important and the shock 
develops, we may regard it as a standing wave of n=l (the funda- 
ment^). Thevelocity of the wave propagation in the medium is 
(see |Painll2005h 

Cw = ojL/tv (56) 

where L is the extent of the domain (0.5 cm) and co the wave fre- 
quency (27r/Ar, AT its period). 

The initial condition for the velocity is to be zero everywhere. 



we therefore choose an arbitrary position in space (x = 0.01 cm) 
and measure the time it takes to Ux to reach its minimum before ris- 
ing again. This corresponds to AT/4. We then compute the value 
of c„. We limit ourselves to this early stages only to avoid non- 
linear effects. The values we measure are 200.774 cm s~^ for the 
50x50 simulation, 200.793 cm s"^ for the 100x100 simulation and 
200.797 cm s~^ for the 200x200 one. These are very close to the 
expected value for a gravity wave with speed given by y/gH — 200 
cm s^^ in shallow water approximation, although of course we are 
not modelling shallow water here but shallow compressible gas. 
However at this modest amplitude the compressibility has only a 
small effect. In the case of the strong shock we find 192.197 cm 
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Figure 5. Relative variation of E for the wave simulation at resolution 
200x200 (section irm . 



Lx, t = 0.00e + 00 s 




t = 7..5"e + 0:; s 




Figure 6. Snapshot of the Rossby adjustment simulation 200x200 (section 
\4.2l at times t = and 1.52 X 10~^ s. Units of the velocity color scale in 
cm s~^. The transient waves are visible. 



s~^, which is less accurate, but the perturbation is higher and the 
shock sets in at the first crossing, therefore the linear approximation 
is definitely not valid. 



4.2 Rossby adjustment problem 

We also simulate the Rossby adjustment problem. This is a 2D 
problem (one vertical and one horizontal dimension) similar to the 
dam break, with the addition of the effects of the rotation of the 
reference frame. The fluid is assumed to be confined and at rest till 



at f = s it is let free to move. The initial conditions correspond 
to a central "bump" in the fluid which will try to spill laterally un- 
der the action of pressure/gravity. Coriolis force will oppose this 
motion and the fluid should adjust to an equilibrium configuration 
with a sloping interface that extends for ~ 6iiR in th e horizon- 
tal (x ) direction; i?R is the Rossby radius defined as (see lPedloskvl 
[l9^: 



Rr = 



(57) 



where SI is the angular velocity of the reference frame; note that 
below we use the Coriolis parameter / — 2Q. 

The set-up for these simulations is: g — 10^ cm s~'^, Pt = 
10* erg cm~^ P. = (e - 1)Pt and PT^ax = 2 x 10*^ erg g-\ 
We add an initial perturbation of the fluid temperature as 

(2r;//fniax - 1) oxp [{x - xo)/5x] 



(58) 



1 + exp [{x — xo)/5x] 

This perturbation goes from to 2ri/Hm-^^ — 1 over an interval 5x, 
being rz/f/max — 0.5 at a::o. a;o is chosen to be at 75 per cent of the 
domain, 5x is 0.2 per cent of it and 277/J/niax = 0.5. We measure 
the average height Ho at xq. With these choices Hrasx = 4 x 10^ 
cm, Hq = 3 X lO'' cm and 77 = 10'^ cm. Again, we simulate only 
half of the domain (480 km, 240 km) to save computational time 
and then mirror the results and we try three resolutions of 50x50, 
100x100 and 200x200 at / = 10"^ s"^ (Fig.Hll. 

4.2.1 Adjustment 

Our simulations never reached the steady state, which was to be ex- 
pected given the reflecting boundary conditions and the fact that the 
time to relax can be extremely long (see Kuo & Polvani 1997). In 
order to have a measure of the asymptotic configuration we average 
the profile of the simulations after t = 5 x 10'^ s when only steady 
gravity waves are left. The theoretical prediction for the asymptotic 
shape of the profile in shallow water approximation should be (see 
IKuo& Polvani 1997; Boss & Thompson 1995'): 

f ^i-\/l^^exp[ {x-x^)/Rx] xi^x, 
H2~ J^Aexp[~{x-x^)/R2] x^x^ 



_j Aexp[ {x — X!i)/Rl] X ^ Xa 

" \ Aexp [-{x - Xs^)/R2] x'^x^ 
=0 



(60) 
(61) 



but note of course that this is not completely applicable here as 
our gas is compressible. H\ = Hq + rj and H2 ~ Ho — 77 are 
the maximal and minimal initial heights, Ri = ^/gHi/ f, R2 — 
\fgH2l ^ , the Rossby radii of the two heights, and x^ = R\ — R2, 

A = fx^. 

As it can be seen from Fig. |7] the approximation of the the- 
oretical results improves quite well with the resolution. This is to 
be expected, since the relevant length scale _Rr corresponds to only 
~ 1.8 grid cells in the 50x50 simulation, while it improves to ~ 3.6 
in the 100x100 one and to 7.2 in the 200x200 one. We measure 
the root mean square differenc^EI between the theoretical predic- 
tion for the profile and the numerical results: it is 1.20 x 10~^, 
3.92 X 10"^ and 2.22 x 10"^: definitely improving. Also the 
approximation of the value of Xa increases: relative accuracy is 



We do not use the relative difference to avoid divergences when the the- 
oretical value is 0. 



© RAS. MNRAS OOO.mfTS] 



10 J. Braithwaite and Y. Cavecchi 



Asymptotic adjustment profile 



0.5 
0.0 
-0.5 
-1.0 










"" ~ — -"'^ 


~^ "~- ^ 




- 
_ 




- 














- 


^ "A 


- 




\ A 








- 




v\ 








- 




























\ 




























1 


1 ... 1 ... 1 . 





1.0 - 



0.5 



-0.5 



-1.0 



-2 2 4 

[RH=1.72e + 06 cm] 



Asymptotic adjustment Uy 




x/R„ 



-2 2 4 

[RB=1.72e + 06 cm] 



Figure 7. Left panel: comparison of numerical results (for the Rossby adjustment simulation (section Pt.2.U at resolutions 50x50 dashed line, 100x100 dotted 
line, and 200x200 solid line) with theoretical prediction (equation [59] dot dashed line) for the height profile. The abscissas are rescaled with respect to the 
and ai'e centred on xq , see eguation lSSl where the initial configuration height was Hq . Right panel: For the same simulations, a comparison of numerical 
results (dashed, dotted and solid lines) with theoretical prediction (equation 1601 dot dashed line) for the Uy profile. 
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Figure 8. Time Fourier analysis of the height at a fixed horizontal position 
(x=2.1 X 10^ cm) for the Rossby adjustment simulation / = 10^'^ Hz at 
resolution 200x200 (section l4.2.1t . The vertical line indicates the frequency 
corresponding to f /2tt. 



63.88 per cent (50x50), 60.30 per cent (100x100), and 53.50 per 
cent (200x200). Fig.|7](right-hand side) confirms this trend: the root 
mean square difference between the theoretical prediction equation 
[eolfor Uy and the numerical results is 6.42 x 10"^ 3.68 x 10"^, 
and 1.99 x 10~^. Note that due to the diffusive high-order nature 
of the code, we cannot reproduce the sharp peak in the theoretical 
prediction and this explains the higher discrepancies. Anyway, the 
approximation of the maximum value improves steadily: relative 
accuracies are 58.48 per cent (50x50), 39.02 per cent (100x100) 
and 17.78 per cent (200x200). 

A Fourier analysis of the height of the surface of the gas shows 
the presence of strong oscillations in addition to red noise at low 
frequencies. The peaks start above f /2tt, indicated as vertical line 
in Fig.jS] which is in good agreement with the theoretical dispersion 



t= 0.0 f"' 



t= 2.3 f"= 



t= 5.3 f"- 



t= 8.4 f"= 



-6 -4 -2 2 4 6 

x/R„ [RB=3.46e + 06 cm] 

Figure 9. Time evolution of the height profile for the Rossby adjustment 
simulation with / = 5 X 10"'* Hz at resolution 200x200 (section l4.2.U . 

relation for the waves: (see |Pedloskvlll987h : 

where A is the wavelength of the wave. 

Finally, as an example, Fig.|9]shows the time evolution of the 
height profile for the simulation at resolution 200x200. 

4.3 Kelvin-Helmholtz instability 

This is a shear instability where two fluids are moving parallel to 
each other with different velocities. We ran a simulation at resolu- 
tion 200x200x16 with a central section (accounting for one third 
of the volume) moving with a velocity Uy — 5 x 10~^ cm s~* 
and the remaining two-thirds of the domain moving with equal and 
opposite velocity. Initially, the temperature and P» are both con- 
stants, and P« is equal to Pt, so that a little under one scale height 
is modelled. We follow the locations of the two fluids with the aid 
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Uy, t = ;).COc + CQ s 



ptracor, t = 8. 866-04 s 







t = B.B6c-G'1 



Figure 10. A simulation of the Kelvin-Helmholtz instability in a square box with resolution 200x200x16 (section [4.31 . In the initial conditions there are 
perturbations at live different wavelengths: A = 1/1, 1/2, 1/3, 1/4, 1/5 cm. Top-left: Initial conditions, Uy. The other three frames are at t = 8.86 X 10~* 
s: top-right: passive tracer used to follow the locations of the two fluids, bottom-left: Uy, and bottom-right: Ux - (Note that these snapshots are taken at a later 
time than the maximum time of Fig. ll2t 



of a passive tracer variable. Finally, to get the instability started we 
give the fluid an initial kick, giving the x-component of the veloc- 
ity with perturbations of the form — 50 cos(27rj//A) cm s^^. In 
the following example, five wavelengths are perturbed at once, the 
largest five wavelengths fitting into the domain. The output of these 
simulations is plotted in Fig.llOl 

In this simulation, mass is conserved to machine precision as 
usual, as in the simulation mentioned previously. As for momen- 
tum conservation, which was not tested previously, the fractional 
change in total momentum in the y-direction is plotted in Fig. [TT] 
It is conserved here to within about one part in 10^. 

We now measure the growth rate of the instability taking the 
Fourier transform of in the y direction along the line x — 1.7 x 
10~^ cm (i.e. at the initial position of the interface between the two 
flows). The initial perturbation should grow at a rate exp (ujt). 

Under the assu mptions of no strat ification and incompressibil- 
ity uj is given by (see lChoudhuri|[l998h : 



100 X Apyp, 




10 20 30 

Time [l.e— 4 s] 



40 



Figure 11. Relative variation of py for the Kelvin-Helmholtz instability 
simulation at resolution 200x200x16 (section|43)- 



2n 

— \lPiP2 



H,2 



Pi + p2 



(63) 



Although these assumptions are not trickily true for our simula- 
tions, still they are good approximations and indeed looking at the 
amplitudes of the five perturbed modes, we see that (as is already 



obvious from Fig.lIOt the shortest wavelength grows the fastest, as 
predicted - see Fig. 1 121 
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Figure 12. Time evolution for the powers of A = 1 / 1 ... 5 cm ^ for the 
Kelvin-Helmholtz instability simulation at resolution 200x200x16 (section 
[431 with seeds at different As (1. solid, 1/2 dotted, 1/3 dashed, 1/4 dash 
dotted and 1/5 dash dot dotted). The smallest A (1/5) is the first to grow. 
(Note that the snapshots of Fig. |10| are taken at a later time) 



4.4 Inverse entropy gradient instability 

A common case in astrophysics is when thermal conduction is 
not enough to bring heat from a lower layer to an upper one. 
This situation leads to an inverse gradient of entropy and conse- 
quently to convection. This kind of instability can be thought of 
as a generalization of the standard text book Rayleigh-Taylor in- 
stability and the driving force is still basically buoyancy, the main 
difference being the compressibility of the gas. The criterion for 
stability in case of adiabatic m otion of the flu id elements is known 
as Schwarzschild criterion (see lClavtonlll984l ): 

ds/dz>0 (64) 

where s is the specific entropy and z the height. In the cr-coordinate 
system this translates in to 

ds/da s; (65) 
In the case of an ideal gas we have 

where 7 = cp/cy, which can also be rewritten, with the use of 
equation [6] as: 

T oc pi-i/^e"/T (67) 



s oc In 



(66) 




2. 5x10"^ 
entropy, t=1.14e + 00 s 



2.5x10"- 
cntropy, t = 2.03e + 01 



Figure 13. Snapshots of the entropy field at t = 0, 1.14, 20.3 s for the in- 
verse entropy gradient simulation 200x800 (section |4!4| . Note that entropy 
is not fully conserved, due to mixing. 



For our simulation, we set 

RT ( P 



1 erg \ 1 erg cm 



1-1/7 



which ensures a gradient for entropy of ds/dcr = O.OO37 > 
in violation of condition ( I65t . Pt ~ 1 erg cm~'^ and P* = 
(e — 1) Pt, so that we simulate 1 scale height. The average sound 
speed is ~ 1.5 cm s^^ . To start the instability we perturb the initial 
velocity field as: 



5 X lO^^^sin f + ) cms ^ (69) 

Z — 1 \ * 
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Figure 14. Simulation of tlie propagation of a plane Alfven wave in the 
veilical direction (section 14.5. U - a time sequence of i?J (left) and 
(right). The first six points in time (solid, dotted, dashed, dot-dashed, dot- 
dot-dot-dashed, long-dashed) near the bottom of the plot are at times t = 0, 
46, 93, 139, 186 and 232 s. The next two - the solid lines half way up, 
are at times t = 487 and 835 s. The final six (with the same line-styles as 
the first six) are at times t = 1171, 1218, 1264, 1310, 1357 and 1403 s. 
Note how the amplitude of the wave increases as lower density is reached, 
as viewed in B* . 



where Ai — 1/i cm and ifi is a set of random phases. The domain 
extent is 1 cm. 

In order to make sure that the motion of fluid elements is as 
adiabatic as possible we include just the hyperdiffusive thermal 
conduction (see sec. 13.4.2V This test is run in 2D with a resolution 
of 200x800 and Fig.lOlshows the initial conditions and the evolu- 
tion of the entropy profile: after ~ 20 s the profile has completely 
overturned. 
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Figure 15. Simulation of the propagation of a plane Alfven wave in the 
vertical direction (section 14.5.1) - the position (in terms of coordinate cr) 
of the peak of the wave, as it propagates back and forth between top and 
bottom. 
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Figure 16. Simulation of the propagation of a plane Alfven wave in the 
vertical direction (section I4.5.U - the speed of propagation of the wave 
(solid fine) compared to the theoretical prediction (dashed fine). The two 
agree very closely, except when the wave bounces off the boundaries at t ~ 
1300 s, when it is impossible to measure the propagation speed properly. 



4.5 Magnetic field tests 

In this section the implementation of magnetic fields is tested, by 
propagating Alfven waves and by modelling the Tayler instability 
in a toroidal field. 

4.5.1 Alfven waves 

We test here the propagation of a plane Alfven wave in the verti- 
cal direction. The initial magnetic field is simply a uniform field 
Bz = Bo, and it is set in motion with an initial velocity field 
Ux ~ uomax(0, (cr — ao)/(l — <7o)), which is just a 'hockey- 
stick' shape with non-zero value at a between ao and 1. We set 
(Jo = 0.92, so we have a kick at the bottom of the domain. Periodic 
boundaries are used in the two horizontal directions; in the vertical, 
we use antisymmetric conditions for and symmetric for B±. 
The computational domain has a height equal to one scale-height, 
and the temperature is uniform. As can be seen in Figs. ll4lto[m 
the wave propagates upwards, growing in amplitude as it does so 



in response to the lower density higher up, reflects from the upper 
boundary and propagates back downwards. We follow the propaga- 
tion through a number of journeys between top and bottom, finding 
that the wave is damped only rather slowly. The speed of propaga- 
tion o f the wave is compared to the local Alfven speed jChoudhuril 
Il998h va = Bo /\/4:TTp in Fig.[T6l where we can see that the agree- 
ment is close. The vertical flux is conserved perfectly in this simple 
setup. 

4.5.2 Tayler instability 

This is an instability of a toroidal magnetic field ( lTavled[l957l 
I1973I) . The free energy source is the field itself and the energy 
is released by an interchange of fluid with weaker toroidal mag- 
netic field B^ with fluid containing stronger toroidal field at greater 
cylindrical radius tu. In a field given in the usual cylindrical coordi- 
nate notation by B^ — Bo'co /'coo we expect the m = 1 azimuthal 
mode to be unstable and the growth rate to be roughly equal to the 
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Figure 17. Simulation of the propagation of a plane Alfven wave in the 
vertical direction (section 14.5.11 - the amplitude of the wave (as measured 
by peak ) against time, showing a gradual decay. This decay is reduced 
at higher resolution and/or smaller diffusion coefficients. 
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Figure 18. Simulation of the Tayler instability (section [4.5.2) . The lines 
show the positions of co-moving fluid surfaces as they intersect the y = 
cm plane at five different times: up^t=-Vi, 1.55, 2.10, 2.63 and 3.20, repre- 
sented respectively by the sohd, dotted, dashed, dot-dashed and dot-dot-dot- 
dashed lines. The horizontal extent of the computational domain in both x 
and y is from —2 to +2 cm; only the central part of the y = cm plane is 
plotted here as nothing is happening towards the edges of the box. 



Alfven frequency given by cja = 'va/'^ = i3o/(^o\/47rp)- The 
instability can be modelled in a square computational box contain- 
ing a magnetic field of the form = Bq {vj / wo) /{ 1 + exp[(ti7 — 
wo) /Aro]}, the latter function simply being a smooth taper so that 
the field goes towards zero at the edge of the box. We use the same 
boundary conditions as in the previous case. 

The horizontal size of the box is 4x4 cm^ and we set vjq = 
4/3 cm, Atu = 0.12 cm, Bo =0.1 G, the temperature at the be- 
ginning is uniform and of value RT = 1 erg g^^ and we set g = 1 
cm s~^, so that the scale height = 1 cm. The vertical extent 
of the model is Q.QlHp, which means that all vertical wavelengths 
are expected to be unstable - a strong stratification stabilizes the 
longer wavelengths. A resolution 72x72x72 is used. The code suc- 
cessfully reproduces the instability at all expected wavelengths, and 
the growth rate measured corresponds to that expected (Fig.|18t. 



5 SUMMARY 

We have described a numerical magnetohydrodynamic scheme de- 
signed to model phenomena in gravitationally stratified fluids. This 
scheme uses the a-coordinate system, a system which basically em- 
ploys pressure as the vertical coordinate. In order to do this the code 
assumes hydrostatic equilibrium in the vertical direction. Our code 
is tailored for problems that fulfil the following conditions. Firstly, 
the fluid under consideration should have strong gravitational strat- 
ification, with much greater length scales in the horizontal than in 
the vertical direction (perhaps greater than the scale height H^). 
Secondly, the timescales of interest should be longer than the ver- 
tical acoustic timescale (Hp/cs). Finally, in the magnetic case, a 
high plasma-/3 is required so that Alfven-wave propagation in the 
vertical direction does not limit the timestep. 

The code has been successfully validated. It is capable of re- 
producing very different phenomena like the Kelvin-Helmholtz and 
the inverse entropy gradient instabilities, waves and shocks, the 
Rossby adjustment problem as well as the propagation of Alfven 
waves and the Tayler instability. The code converges well to the 
analytic solutions and conserves mass, energy and momentum very 
accurately. 

This demonstrates our numerical scheme to be both highly 
flexible and the natural choice for many astrophysical contexts such 
as planetary atmospheres, stellar radiative zones, as well as neu- 
tron star atmospheres. We have already used it to investigate flame 
propagation in Type-I X-ray bursts on neutron stars: results from 
this study will be presented elsewhere (Cavecchi et al. in prep.). 
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